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Abstract 

The precise location of the water ice condensation front (snow Hne) in the protosolar 
nebula has been a debate for a long time. Its importance stems from the expected sub- 
stantial jump in the abundance of solids beyond the snow line, which is conducive to 
planet formation, and from the higher stickiness in collisions of ice-coated dust grains, 
which may help the process of coagulation of dust and the formation of planetesimals. 
In an optically thin nebula, the location of the snow line is easily calculated to be 
around 3 AU, subject to brightness variations of the young Sun. However, in its first 5 
to 10 million years, the solar nebula was optically thick, implying a smaller snowline 
radius due to shielding from direct sunlight, but also a larger radius because of viscous 
heating. Several models have attempted to treat these opposing effects. However, until 
recently treatments beyond an approximate 1-i-lD radiative transfer were unfeasible. 
We revisit the problem with a fully self-consistent 3D treatment in an axisymmetric 
disk model, including a density-dependent treatment of the dust and ice sublimation. 
We find that the location of the snow line is very sensitive to the opacities of the dust 
grains and the mass accretion rate of the disk. We show that previous approximate 
treatments are quite efficient at determining the location of the snow line if the energy 
budget is locally dominated by viscous accretion. Using this result we derive an an- 
alytic estimate of the location of the snow line that compares very well with results 
from this and previous studies. Using solar abundances of the elements we compute 
the abundance of dust and ice and find that the expected jump in solid surface density at 
the snow line is smaller than previously assumed. We further show that in the inner few 
AU the refractory species are also partly evaporated, leading to a significantly smaller 
solid state surface density in the regions where the rocky planets were formed. 
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1. Introduction 



The efficiency of the formation of rocky planets and the cores of gas giant planets 
depends sensitively on the amount of solids present in each location of the solar nebula. 
In most models of planet formation the solar nebula has been assumed to resemble the 
structures proposed by | Weidenschilhng] ( | 1 977| l and |Hayashi| ( |198l"| ). In this model, and 
subsequent incarnations of it, the surface density of the gas is approximately given by 

2f,rV)- 1700(^) gram/cm^ (1) 
(see Thommes and Duncan 20061) where the acronym MMSN stands for "minimum 



mass solar nebula", a name by which this model is often identified. In the inner regions 
(small values of R) of the nebula the dust grains are refractory in nature, while beyond 
some radius Ri^e water ice is present, most likely in the form of ice mantels surrounding 
the refractory grains. In the Hayashi MMSN model this location is assumed to be 
Rice - 2.7 AU, which is consistent for temperatures expected in an optically thin nebula. 
Often this water ice sublimation boundary radius is called the 'snow line' . Since ices 
may contain a considerable mass compared to the mass in refractory grains, the total 
surface density was assumed to jump atR- Rice- 

EMmsN(^) = 7.1Fi,e ( — ) gram/cm^ (2) 

with 

1 , for < R.ne 



(3) 

4.2, for R > Rice 

( |Thommes and Duncan[ [2006) 1 

This model serves in many planet formation models to compute the midplane gas 
and solid density and the isolation mass of planetary embryos. However, it has been 
shown that the outcomes of planet formation models depends critically on the value 



of Rice and the strength of the ice jump Fi^e (see e.g. Mordasini et al. 2009 1. Also 



dust coagulation and planet formation and migration models are very sensitive to these 



parameters (e.g. Brauer et al. 2008 Ida and Lin 2008 1. Simply taking the values to be 



Rice - 2.7 AU and Fice - 4.2 (for R > Rice) is probably too simplistic and may lead to 
major errors in predictions of the progress of the planet formation process. 

Besides the amount of solid mass available also the strength of the bonds between 
dust grains is of importance to their aggregation behavior The presence of ice coat- 
ings on grains can significantly increase their 'stickiness', allowing grains to coagulate 
more easily. Perhaps even more importantly the bonds formed are stronger preventing 
aggregates that formed to be destroyed again by collisions. 

Over the last decade or so, there has been enormous progress in astronomy in the 
understanding of the structure and evolution of 'protoplanetary disks', i.e. the dust+gas 
disks surrounding many very young stars. These disks very likely give a good impres- 
sion of what our own solar nebula looked like 4.567 billion years ago, and the knowl- 
edge gained in that field can help getting a better handle on the structure of the solar 
nebula and the location of the ice condensation front. Indeed, modeling tools that were 
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originally meant for modeling protoplanetary disks and comparing model predictions 
to observations have been used to make models of the solar nebula and the snow line. 



For instance, the 1 + lD disk structure model of Hure (20001, in which the equations 



are first solved vertically (ID) and then connected radially (1 + lD), has been used by 



Hersant et al. (2001jl to study water in the solar nebula, in particular concerning the 



issue of the D/H ratio in meteorites. |Davis| ( [2005b| l made a model of the surface den- 
sity of the solar nebula disk as a function of radial coordinate and accretion rate. He 
used this model to study the ice condensation front ( |Davis[ |2005a) . This model in- 
cludes both the effect of irradiation of the disk by the young sun and the heating of the 
disk near the midplane through viscous dissipation of potential energy (i.e. accretional 
heating). Like in the Hersant study, these models are 1h-1D models, i.e. they are 1-D 
vertical models of density and temperature as a function of z. The Davis models in- 
clude radiative transfer using a variable eddington factor method. [^ Lecar et al. (2 006) l 
also modeled this, using an updated version of the |Chiang and Goldreich,(1997^ two- 
layer disk model, i.e. slightly simpler than the 1h-1D modeling of Davis. Recently, 



Dodson-Robinson et al. p009] l combined complex models of ice formation and vis- 
cous evolution of the disk with a simple radiative transfer approximation to compute 
the evolution of the solid surface density. 

While these modeling efforts have improved the understanding of the distribution of 
ices in the solar nebula substantially, they still rely on rather dramatic simplifications of 
the treatment of radiative transfer The reason for this is very understandable: it is quite 
a numerical challenge to model multi-dimensional radiative transfer in protoplanetary 
disks that are very optically thick and actively accreting (and hence producing heat near 
the midplane). However, in recent years multi-dimensional radiative transfer tools, 
in particular those based on the Monte Carlo method, have improved dramatically in 
speed. Extreme optical depths are now not necessarily a problem anymore. We have 
developed a code, MCMax, that can now handle extreme optical depths efficiently, 
even if many of the photon packages originate from the optically thickest regions of 



the disk (Mm et al. 2009 1. In addition, we have adjusted the code to fully take into 



account the density dependend sublimation state of both the refractory as well as the 



icy material (Kama et al. 2009jl. It is therefore a right time to revisit the problem of 



the snow line with the full force of multi-dimensional radiative transfer and investigate 
(a) where is the ice line as a function of time and (b) what is the gas and solid surface 
density distribution in the disk as a function of time. This is the goal of this paper. 

2. Model 

Our model combines the radiative transfer equipment and basic disk model setup 
described injMmeFaL ( 2009| l. We add here: a self-consistent and automatic computa- 



tion of the radial structure, viscous heating of the disk by accretion and a self-consistent 
treatment of water ice and dust evaporation. 

Right at this point we already note that we will replace 'time' t with 'accretion 
rate' M and study only the inner 20 AU of the disk. The reason is that the long term 
evolution of protoplanetary disks depends very strongly on processes such as photoe- 
vaporation of the outer disk regions (see e.g. Gorti and Hollenbach 2009| l. The inner 



disk regions, say inward of 20 AU, are regulated entirely by the 'feeding' of matter by 
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the outer disk. How this feeding depends on time requires a detailed study of the com- 
plex processes happening at large radii. But given a recipe of feeding, i.e. an accretion 
rate as a function of time M{t), the inner disk can be assumed to be entirely determined. 
In particular if we avoid issues such as ionization degree and possible instabilities in 
the disk (see e.g. Hartmann et al. 2006') and assume that no mass pile-up can happen 
anywhere in the inner disk (this assumption must be relaxed in future work), then the 
stationary disk structure inward of 20 AU is entirely determined by the stellar parame- 
ters and the accretion rate M. For any given star our model thus becomes a 1 -parameter 
model. So instead of having t as our parameter, we will have M as our parameter Gen- 
erally, a high mass accretion rate corresponds to early times in the evolution of the solar 
nebula, while low mass accretion rates correspond to later times. Alternatively, in the 
episodic accretion scenario, low mass accretion rates represent the typical state, which 
is punctuated by relatively short episodes with high mass accretion rates. 

The main improvement of our study over earlier work is that we treat the flow of 
radiation through the disk in a fully self-consistent manner, assuming that dust opacities 
are everywhere dominant over the gas opacities. Our model is axially symmetric, i.e. 
all variables depend on radial coordinate r and vertical coordinate z, but do not depend 
on longitudinal coordinate (p. The movement of photons is, however, fully 3-D, i.e. a 
photon can move also outside of the r, z-plane and has all the 3-D freedom of motion. 
In our Monte Carlo scheme (based on Bjorkman and Wood 2001] ) we split the total 
input luminosity L into 'photon packages', each carrying a luminosity Lpackage - 
L/N. The total luminosity includes the stellar luminosity Lt and the accretional heating 
luminosity Lvisc, i-e. L - L„ + Lvisc- As these photon packages travel through the disk 
they update the local temperature of the dust according to the scheme by |Bjorkman| 
and Wood] ( 2001 | l. Since our model is axisymmetric, these updates have to only be 



done in 2-D, i.e. the temperature only depends on (and will only be updated in) r, z, no 
matter what is the longitudinal location (p of the photon package. So while the radiative 
transfer is 3-D, the resulting structure is 2-D. 

The main difference between models doing radiative transfer in 1-i-lD and our 
model doing full 3D is in the radial energy diffusion. In the 1-i-lD model the energy 
from the star is intercepted by the disk and further diffusion is assumed to be vertical. 
This approximation is expected to be accurate for the regions of the disk close to the 
midplane where the very high densities do not allow much radial diffusion. However, 
radial energy diffusion is important in the surface layers and, more importantly, in the 
innermost regions of the disk, directly illuminated by the central star. Our model natu- 
rally takes this into account and allows us to see in which parts of the disk the vertical 
diffusion approximation is valid. Also, 3D radiative transfer is the only way to properly 



treat the complex shape of the dust condensation front close to the star (see alsO |Kama 
[etiLilSOOQl l. 

For the central star in our model we take the Sun, i.e. the mass = Mq, the 
luminosity L* - Lq, and the stellar radius - Rq. The spectrum of the Sun is 
approximated by a blackbody of 5777 K. Although the parameters of the early Sun 
were likely different from this, we do not expect much difference on the results in this 
paper. 
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2.1. Treatment of viscous accretion 

We assume that viscous heating is described by the so-called c-disk model ( ShakuraJ 
|and Sunyaev||197T) ). In this model it is assumed that the energy produced locally by 
viscous heating is proportional to the local gas pressure through the proportionahty 
constant a. The total energy locally created per unit volume is then given by 

Tviscous = '^-aP{z)mR), (4) 

where P{z) is the gas pressure, and Q.{R) is the Keplerian angular velocity at radius R 
from the central star. If we define a uniform accretion rate M throughout the disk we 
have for the total energy released per unit surface area of the disk per unit time at a 
radius R from the Sun is 



—R-^ 

An 



(5) 



where G is the gravitational contstant, and Mi, and 7?* are the mass and radius of the 
central star. Combining Eqs.|4]and|5]we can solve for the surface density of the disk 
for any given value of a. In this paper we will consider two values of a, 0.1 and 0.01. 
From Eq. (j5]l it is clear that the total amount of energy released is independent of the 
structure of the disk. This speeds up convergence since there is no feedback from the 
disk structure on the total energy that is released. 

We basically treat the energy released by viscous stress in the same way as the 
energy released by the central star Thus a fraction of the photon packages released in 
the Monte Carlo radiative transfer process is now released from the inner regions of 
the disk. Where previous studies that take into account radiation from viscous stress 
rely on assumptions on the radiative transfer in the optically thick inner regions, we 
treat the radiative transfer fully consistent in the framework of Monte Carlo radiative 
transfer as outlined byjMin et al.|(|2009!l. The total accretional heating luminosity of the 
disk is given by the integral of Eq. (j5]l from the stellar surface, to the outer radius, 
.Rout- Assuming Rq^ » R* this leads to 

_ 1 GM.M 

i-visc - X ^ ■ W 

Z Ri, 

Note that a large fraction of this energy is most likely released inside the dust sublima- 
tion radius and thus can easily escape the disk. This energy does not influence the dust 
temperature and is thus not taken into account in our computation. In our method the 
energy produced through viscous heating is directly deposited from the gas onto the 
dust grains. This assumes a strong coupling between gas and dust, which is an appro- 
priate assumption in the regions where the densities of both gas and dust are high. In 
the regions where a large fraction of the dust is evaporated due to high temperatures, 
we limit the amount of energy deposited on the grains by viscous heating by taking 
into account the number of collisions between the gas and the dust and assuming a gas 
molecule can at maximum transfer its total kinetic energy onto the grain in such a col- 
lision. This limits the total energy that can be transferred to the dust per unit time. We 
find that indeed in the optically thick regions this limit exceeds the available energy 
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Table 1 : Parameters of the three different dust models considered. The composition is given as mass frac- 
tions. The gas to solid ratio is the ratio between the mass in gas and the total available mass to form solids. 
Ficc is the fractional increase in the mass available to form solids when water ice can condense. The value of 
the Rosseland mean opacity, a:^:, is taken at 160K is computed considering all dust species including water 
ice. 



by many orders of magnitude. However, in the inner regions, this limit prevents the 
small amount of dust grains available to be overloaded with accretional energy from 
the dense gas surrounding it. The dust grains are assumed to radiate the energy they 
absorbed with a wavelength distribution according to the local dust temperature after 
which the regular procedure of Monte Carlo radiative transfer is proceeded. 

In the dusty regions of the disk the optical depth from the midplane, where most 
energy is released, to the outside is extremely large. Therefore, a photon package un- 
dergoes a large number of interaction steps before it leaves the system (on the order 
of several million interactions). If all these interaction steps would be computed sep- 
arately, the computation time required would be on the order of a minute per photon 
package, putting strong constraints on the number of photon packages that can be used. 



Min et al. ( 2009| introduced an analytical solution for the random walk regime of the 



radiative transfer that can be applied in the optically thick regions and reduces compu- 
tation times by orders of magnitude. This algorithm therefore allows the use of a large 
number of photon packages and thus to obtain proper statistics everywhere. 

2.2. The dust and ice composition and opacities 

In order to estimate the abundances of the various dust components we use here the 
Solar abundances of the elements as determined by Grevesse and Sauval ( 1998| l. The 
dust mixture is then constructed as follows: 

• Since iron sulfide is the dominant form of (solid) sulfur in meteorites in the solar 
system it is first assumed that all the available sulfur is in the form of FeS. This 
takes approximately half the available number of Fe-atoms. 

• Then we assume all Mg and Si go into silicates together with the remaining Fe- 
atoms. This creates silicates with an average composition of roughly MgFeo sSiOs.s. 
We assume this to be in the form of amorphous silicates. This takes approxi- 
mately 15% of the available O-atoms. 
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• Next we put a fraction w of the available C-atoms in the form of carbonaceous 
dust grains. The remaining C-atoms are assumed to form CO. This takes away 
up to ~ 50% of the available number of 0-atoms depending on the value of w. 

• The remaining O-atoms are assumed to form water ice everywhere in the disk 
where the temperature allows it. 

From the above procedure we also can compute the gas to dust ratio. We find that the 
resulting gas-to-dust ration can deviate up to a factor of 2 from the canonical value 100. 
The question remains what the most realistic value of the carbon abundance in dust, w, 
is. We take as a standard the value w - 0.5, based on in-situ measurements of the dust 
in comet Halley (Geiss 1987| l. However, the carbon fraction in dust could have been 



higher in the early Solar nebula since the interstellar solid carbon abundance is much 
higher. Also, in the inner regions, at temperatures above ~ 1000 K, carbon grains will 
be destroyed by combustion so lower values of w are to be expected. The dust/ice/gas 
mixture resulting from this procedure for different choices of w are summarized in 
Table [U 

In order to compute the opacities of the dust species as a function of wavelength 
we need to know the refractive indices of the grains. For the silicates we use the mea- 
surements by Dorschner et al. ( 1995 1 and Henning and Stognienko ( 1996| l, for the iron 



sulfide we take the laboratory data from'Begemann et al.'('1994), for the carboneceous 
material we take the measurements from Preibisch et al._(1993) , and for the water ice 
we take the data from |Warren| ( | 1984 1. In order to convert the refractive index data to 



opacities we use the method by 'Min et al. ( 2005 | l to simulate irregularly shaped parti 



cles. The particle sizes are assumed to follow the so-called MRN distribution ( ,Mathis 



et al. \911) . We thus assume that no significant grain growth has occurred with respect 



to the grains in the interstellar medium. We will come back to this assumption in the 
discussion of the results. 

2.3. Treatment of evaporation and condensation 

The temperature at which a species evaporates depends on the vapor pressure of 
the material. We treat the evaporation and condensation of the dust species following 



Kama et al. ( 2009| l. In order to converge the structure of the inner regions pCama et aL 
( 20091) had to put constraints on the vertical and radial gradients of the fraction of the 



material in the gas phase. This is to avoid instabilities caused by the sudden jump in 
opacity at the inner edge of the disk. Since the ice is located in the inner regions of the 
disk, the ill conditioned convergence is no issue here. Therefore, for the water ice we 
do not put any constraints on the gradients of the condensed fraction. It turns out that 



convergence is reached relatively fast and no instabilities like those reported by Kama 



et al. ( 2009 1 are found for the ice sublimation zone. Although the density structure 
does not change significantly anymore after ~ 30 iterations, we ran all models up to 
100 iterations to make sure also the slowly changing structures are converged. 

To simplify the evaporation of the refractory dust species we take the evaporation 
parameters for amorphous olivine for the entire dust mixture of silicates, iron sulfide 
and carboneceous material. For the iron sulfide it might be argued that this is not a 
correct assumption since the evaporation temperature of FeS is around 680 K. However, 
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Figure 1: The areas in tlie disk where water ice is present for various mass accretion rates. The case in 
which half of the available C atoms are in the form of solid carbon is presented for a value of the viscosity 
parameter a = 0.01. Increasing intensity corresponds to the region where up to 25%, 50%, 75% and 100% 
of the available ice is condensed. The white line corresponds to the location where the radiation from the 
star hits an optical depth of one. 



we assume that above this temperature the iron condenses as metallic iron, contributing 
similarly to the opacity of the grains. The evaporation characteristics of both ice and 
oUvine are based on sublimation data from |Pollack et al.] ( |1994[ l. 

2.4. Convection 

In the implementation of radiative transfer and the computation of heat balance in 
the disk we employ, the effect of convection is not taken into account. However, we 
will show that convection is an important way of transporting energy generated in the 
midplane of the disk by accretion to the surface. Therefore, we have implemented a 
rough treatment of convection in the radiative transfer code to test the effects on our 
results. In the computations below, convection is not taken into account unless stated 
otherwise. 

Convection becomes an important way of transporting heat when the radiative tem- 
perature gradient becomes too large. According to the Schwarzschild criterion for 
convection this happens when 

dinP 

where T and P are the gas density and pressure, and is the adiabatic limit which 
we take to be that of a diatomic gas: V^d - 2/7. Our rough implementation to correct 
the obtained temperature structure for convection is simply as follows. We start at the 
surface of the disk, where the temperature gradient is very small and work our way 
to the midplane. At each location the criterion according to Eq. |7] is determined. If 
the temperature gradient is too steep the temperature gradient is adjusted such that 
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Figure 2: Same as Fig.^but for the case where all carbon is in the form of solid carbon. 



Eq. |7] becomes an equality. In some cases this causes the temperature structure close 
to midplane to be adjusted such that the midplane temperature is reduced. This way 
we can model the effects of cooling of the midplane by convection although we should 
stay aware of the approximate nature of the implementation. 



Our implementation of convection is very similar to the one by D'Alessio et al. 
( |1998[ ) with the difference that we do not consider the turbulent flux, and we take 
the convection efficiency, their parameter f , to be unity. In the relatively less dense 
regions of the disk transfer of energy through radiation dominates and the criterion 
Eq. (jTjl is fulfilled. In this case we do not alter the temperature structure. However, in 
those regions of the disk where radiative transfer cannot fulfill Eq. (j?]) we do adjust the 
temperature structure to be adiabatic. This means that convection becomes important 
only in the more massive disks where radiative energy transfer is difficult. 



3. Results 

3.1. Location of the snow line as a function of M 

In Fig. [T]we display the locations in the disk where water ice can exist for various 
values of the mass accretion rate. It is clear that the location of the snow-line depends 
strongly on the value of M. This is caused by two main effects. The first is that when 
M is increased, the energy released by viscous stress is increased. Secondly, a larger 
value of M causes the disk to be more massive (see also Fig.|4]i and thus the greenhouse 
effect, keeping the energy at the midplane, is much stronger (as explained below, see 
also Eq. ([8])) and the midplane temperature is increased further 

Going vertically up in the disk, the effect from viscous heating becomes less im- 
portant and the snow line moves in significantly. This continues up to the point where 
radiation from the Sun becomes the dominant energy source and causes the snow line 
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Figure 3: Same as Fig.^but for the case where all C atoms are in the form of CO and a = 0.1. 



to move out again. The fact that the location of the snow line in the upper, optically 
thin region of the disk is not constant is caused by the pressure dependent evaporation 
temperature of the ice in combination with a decreasing density when going higher up 



in the disk. This shape of the ice region is consistent with the findings from Davis 



( 2005a I but the snow line in our computations lies significantly further out. The rea- 
son for this is most likely predominantly caused by the different opacity of the grains 
which is assum ed. [ Davis ( 2005 a[ ) use for the opacity the analytical expression from 
Ruden and Lin| ( |1986| l. As outlined below the most important parameter determining 
the midplane temperature is the Rosseland mean opacity at 160 K divided by the gas to 
solid ratio, since this determines how efficiently the energy created by viscous stress is 
kept in the midplane at the location of the ice evaporation front. In the study by [Davis 
( 2005 a| this parameter is krIJ - 2.7cm^/g, while in our study, taking a realistic dust 
grain model, we find that lies in the range 5.9 to 11.1 cm^/g (see Table[T]). Thus in 
our model the energy generated by viscous heating is more efficiently captured in the 
midplane regions. 

In Figs. [2l and [3] we show the extreme cases of the range of dust parameters we 
studied. Fig.|2]shows the case where all carbon atoms are locked into solid carbon, i.e. 
the opacity is at its maximum. In this case the greenhouse effect, heating the midplane 
efficiently is at its maximum. The other extreme, presented in Fig. [3] shows the case 
where all carbon is in CO, i.e. ffie opacity is minimum, and the turbulence is very high, 
a = 0. 1 . This high value of a implies that one can get the same mass accretion rate with 
a smaller surface density. Thus for a given value of M the surface density drops and 
ffie greenhouse effect is smaller In this case the energy can easily escape the midplane 
and the midplane temperature is low leading to a snow line close to the star 
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3.2. Analytical estimate for the location of the snow line in the midplane 

We take the analytic estimate of the midplane temperature from Hubeny ( 1990| l. 
This only includes viscous heating in the midplane and results in a midplane tempera- 
ture, assuming high optical depths, of 



1287To-R^f ' 



(8) 



where kr is the Rosseland mean opacity of the gas/ice mixture and / is the gas to solid 
ratio such that the total optical depth through the disk is equal to Ttot = ArsEgas// = 
fs^soiid- The values of kr and / for the models we computed are listed in Table[T] The 
surface density of the gas, Sga,, can be expressed in terms of T and a using Eqs. Q 
and dSll. The integral of Eq. (j4| in vertical direction is equal to the total energy given in 
Eq. (pji. Assuming R » and taking the temperature to be the midplane temperature 
at each height in the disk we get (see also |Dullemond et"ar 2007 1 



finipM GM. 



3nakhT 



(9) 



where fi is the average molecular mass and nip is the proton mass. Substituting Eq. (j9^ 
into Eq. ^ we get 

rGM* 



5 ifinipM^KR 



USn^akbO-f 



Ri 



(10) 



Taking the above we find that the snow line is located at 



C 



(11) 



where kr is the Rosseland mean opacity of the gas/ice mixture computed at 160 K (see 
Table [TJ, and C is given by 



3/2 



c = 



3yump (GM*) 

ice 



= 1.7- 10 



22^j^5/2 



(12) 



We find that this formula agrees very well with our findings for the range of parameters 
considered here (within a few percent). The formula becomes less accurate when the 
snow line is close to the star, and clearly the equation fails when the radiation of the 
Sun is the dominant energy source in the midplane. However, this is never the case in 
the models considered here. For easy computation, the value of C = 3.5 ■ 10'"* can be 
used in Eq.[TT]when M is given in Mg/yr and kr is given in cm^/g to give the location 
of the snow line /?ice in AU. 

Using Eq. [TT|provides a way to scale our results to those obtained by Davis ( 2005 a| l 
taking into account the correct values of kr, /, and M. Doing this we find that the 
agreement is excellent, and thus the differences are indeed due to the diff'erent values 
of krI f. We also compared our results to those computed by Garaud and Lin ( 2007| l. 
The values for the location of the snowline as computed using Eq. (Ill are compared to 
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This study (w = 0.5, a = 0.01) 
1-10-^ 770 
1-10-^ 770 
1 ■ lO-*^ 770 
1 ■ 10-'^ 770 



16.1 
5.8 
2.1 
0.7 



' 16 

- 6 
2.2 
0.8 



This study (w = 0, c = 0.1) 
1 ■ 10-^ 59 
1 ■ 10-^ 59 
1 ■ 10-** 59 
1 ■ 10'^ 59 



9.1 
3.3 
1.2 
0.4 



- 9 
3.2 
1.2 
0.5 



Table 2: A comparison of the location of the snow line as obtained using Eq. (TT 



computations. The values from the studies by Davis (2005a I and |Garaud and Lm 



and using numerical 
j2007 t are read from 



the figures. The value of KuKfa) for the study by Davis (2005a I is an estimate since they use a dift'erent 
description for the viscosity which corresponds approximately to ours with a = 0.01. For the study by 
[Garaud and Lin (^2007; we used their Eq. (2) to determine the value of at 1 60 K. The values of Rice for the 
numerical computations of this study correspond to the location where 50% of the ice is condensed. 



those obtained by Davis ( 2005a') and 'Garaud and Lin ( 2007 1 in table |2] Note that Davis 
( 2005 a| l adopted a different description for the viscosity of the disk, which explains 
partly the differences. The overall agreement with the approximated and numerically 
computed values implies that estimates of the location of the snow line using approxi- 
mate methods are in general correct when appropriate opacity values are used. 

3.3. Solid surface density distribution as a function of M 

In Fig.|4]we show the surface density as a function of radius for different values of 
the mass accretion rate. Both the total surface density (top panel) as well as the surface 
density in solids (lower panel) are shown. The sudden drop in the total surface density 
at the location where the solids start to become important is an effect of the sudden 
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increase in opacity, causing a steep temperature gradient which translates into a steep 
surface density gradient through Eqs.[4]and[5] The surface density in solids is adjusted 
iteratively until the temperature nowhere in the disk exceeds the sublimation temper- 
ature. This means in the inner regions solid matter is removed until the greenhouse 
effect, keeping the energy locked in the midplane regions of the disk, drops to a level 
where the temperature is low enough. For very high accretion rates, a lot of matter 
has to be removed to let the large amount of energy produced in the midplane escape. 
This balance between condensation and evaporation results in a thermostat keeping the 
midplane temperature at the silicate evaporation temperature over a significant part of 
the inner disk. When the accretion rate is reduced, more matter can condense. How- 
ever, at the same time reducing the mass accretion rate reduces the total surface density 
and thus the amount of matter that can condense. Both these effects play a role in 
determining the surface density in the solid phase. 

At the radius of water ice sublimation we have a gradual increase in the surface 
density in the solid phase. We find that the increase in the solid phase due to water ice 
is Fics = 1.84 for the case where half of the carbon atoms is in the solid phase. This 
increase is quite small compared to the Fi^e - 4.2 taken currently in the literature. We 
have to take into account here that we have only considered water ice condensation. If 
we take in addition CH4 and NH3 ice, which are the most important ice species after 
water ice, we gain about a factor 1.5 (see e.g. Dodson-Robinson et al. 2009 | l so we get 



Fice ~ 2.8. Note however, that these ices condense at lower temperatures, so somewhat 
further out then the water-ice condensation radius we compute here. 

The surface density is in principle a function of the mass accretion rate, the value of 
a and the opacity of the grains. However, in the innermost region the dust evaporates to 
reduce the greenhouse effect such that the midplane temperature drops below the dust 
sublimation temperature. In this regime the solid surface density is adjusted to reach 
a given vertical optical depth, krZ^^sI f, roughly following Eq. ([8]). The maximum 
vertical optical depth that can exist before the midplane gets too hot is thus only a 
function of the total energy released in the midplane, and thereby only of M (see Eq.pb. 
The amount of mass that can exist in this region is robustly determined through the 
mass accretion rate and the opacity of the grains and is independent of the gas density 
or the description of the accretion mechanism. Therefore, in this region we can easily 
scale the solid surface density as a function of M we find in this study to different 
values of the opacity of the grains when, for example, they grow to larger sizes. 

3.4. Effects of convection 

When studying the temperature structure in the disk in detail we find that often the 
temperature gradient is steeper than the adiabatic temperature gradient (see Eq|7]l. In 
these cases the midplane will cool by convection. If we correct for the convective cool- 



ing in the way discussed in section 2.4 the midplane temperature decreases in some 
cases significantly. We find that this only plays a role when the midplane temperature 
is dominated by accretion and in addition, is most important when the temperature is 
above ~ 500 K. Thus, the location of the snow line is only mildly affected. For the case 
of M = 10"^ M^/yi we find that due to the decreased midplane temperature the fraction 
of condensed solids can increase locally by a factor of 5 with respect to the case where 
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Figure 4: The surface density as a function of radius for different values of tlie mass accretion rate. The 
top panel shows the total surface density, while the lower panel shows the suii'ace density in solids. Blue 
solid lines are for M = 10"*Mo/yr, red dotted lines for M = lO^'Mo/yr, green dashed lines are for 
M = 10"** Mo/yr, and purple dashed-dotted lines are for M = 10"' Mo/yr. In all cases half of the available 
carbon atoms are in the solid phase and a = 0.01. 



convection is ignored. For low accretion rates convection becomes unimportant at al- 
most all radii and the disk cools through radiation (see also D'Alessio et al. 1998| l. The 
reason for the dependence of the importance of convection on the surface density of the 
disk is simply due to the fact that less massive disks can easily cool through radiation, 
so convection is not needed to lose the energy from the midplane. In the more massive 
disks it is very difficult for the radiation to escape from the midplane and convective 
cooling is much more efficient. Our treatment of convection is somewhat simplified 
and more advanced methods, also including turbulent fluxes, should be implemented 
in the future to get a clear picture of the exact effects. 

In Fig. [5] we plot the surface density in the solid phase for computations with and 
without convective cooling of the midplane. It is clear that for the high accretion rates 
convective cooling is important to gain extra solids in the region where the terrestrial 
planets are (i.e. in the inner few AU). When convective cooling is ignored we find that 
there is never enough mass in the solids around 1 AU to create the Earth. However, with 
convective cooling we find that there is much more mass. Note that the amount of mass 
in the solids can be increased even further for the high accretion phases by growing the 
grains to somewhat larger sizes, thereby decreasing the opacity and thus the midplane 
temperature. We conclude that convective cooling is an important mechanism to take 
into account when modeling the solid density distribution in the early Solar nebula. 
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Figure 5: The surface density as a function of radius for different values of the mass accretion rate. The solid 
lines are with convective cooling included, the dotted lines are when only radiative cooling is used (same as 
in Fig.|4j. In all cases half of the available carbon atoms can be in the solid phase and a = 0.01. In all panels 
the dashed lines gives the MMSN according to Eq.jljwith and without ice. 



4. Discussion 



4.7. Implications for planet formation 

In this section we compare our results with constraints from the minimum mass 
solar nebula (MMSN). We will use here the 'classical' MMSN as mentioned above. 
There is still quite some debate on what the minimum mass is needed to form our 



solar system (see e.g. Davis 2005b Desch 2007 Crida 2009 1, and most of these 



models find the need for a more massive solar nebula than the classical one we use 
here. Therefore, our choice of the MMSN is conservative, and the computations could 
be repeated to find stronger constraints by using other models. In addition, planet 
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formation might be more dynamic, i.e. planets could accrete mass from different parts 
of the disk through dynamic interactions. 

According to the MMSN (Eq.|2]l in order to form the Earth one needs a solid surface 
density of 7.1 g/cm^ at 1 AU. As can be seen directly from Fig.jsjthis surface density in 
solids is never achieved in our standard model. This is because for the high accretion 
models the midplane gets too hot, so solids have to evaporate to let the energy escape, 
while for the low accretion models there simply is not enough mass at all. However, as 
discussed above, if we take into account convective cooling of the midplane the amount 
of condensed solids can increase by a factor of 5. This way the surface density at 1 AU 
is only slightly below the requirement from the MMSN. Inside 1 AU the surface density 
in solids is never even close to the MMSN, which might pose serious problems for the 
formation scenario of Venus. 

Jupiter is generally assumed to form outside the snow line. Thus, when Jupiter 
formed at its current position, we have as an additional constraint that Ri^e < 5.2 AU. 
For this the mass accretion rate needs to be smaller than M < 8 ■ 10"** Mg/yr. Since 
the midplane temperature at the snow line should be smaller than 160K convective 
cooling is not important, so we can use our standard model. At this mass accretion rate 
the solid surface density is approximately Sjoiid ~ 2 g/cm^, which is slightly lower than 
the mass needed according to Eq. |2] However, typical giant planet formation models 
require up to 3 to 4 times the value of the MMSN (see e.g |Pollack et al. 1996| l. 



We thus find that our computations put constraints on the scenario of the formation 
of the Solar system in the context of the MMSN. However, as we already noted above, 
the MMSN might be an optimistic estimate, and in reality a much larger mass could be 
required. A possible solution to still meet the requirements of planet formation in that 
case would be to lower the opacity of the grains significantly by, for example, growing 
them to larger sizes. 

4.2. Dust properties in the hot midplane regions 

We find that at high mass accretion rates there is a significant region in the disk 
where the midplane temperature is above 1000 K. This is out to ~1.7AU in the case 
where M - 10"^ MJyr, and out to ~0.6 AU in the case of M = 10"^ M^/yv. This is in- 
cluding convective cooling. Thus we find that a considerable region in the early Solar 
system had conditions where thermal processing of the solid state material was very 
efficient. Also, at these temperatures grain coagulation is influenced by sintering (see 
e.g. |Poppe]|2003| l. Sintering is the process where the bonds between grains are molten 
together by the impact of grain collisions or by a long exposure to high temperatures. 
The structures formed under these conditions can be extremely open aggregates when 
the strong bonds created do not allow restructuring or compaction of the aggregates. 
Poppe ( 2003 [ l compute that for aggregates composed of silica spheres exposure to tem- 



peratures around ~ 1250 K is needed for efficient sintering. This is comparable to the 
crystallization temperature of silica found by fabian et al. (20001. Since both sintering 



and crystallization are processes that are caused by the atoms in the lattice wanting 
to move to the minimum energy state, it is not unlikely that the temperatures of these 
processes are similar In that case sintering of Fe/Mg silicates will happen at slightly 
lower temperatures, i.e. around ~ 1000 K. The collision speeds in the midplane regions 
of the disk are too low to cause collisional sintering. 
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As mentioned before, grain growth might help to increasing the solid surface den- 
sity by reducing the opacity of the particles and in this way decreasing the temperature 
in the disk. A low opacity reduces the greenhouse effect and thus allows more grains to 
condense before the midplane temperature gets too high. In fact, when the Rosseland 
mean opacity computed at the dust evaporation temperature of 1500 K reduces by a 
given factor, the maximum allowed surface density increases by this factor. In order to 
test to what sizes we have to grow the grains in order to significantly increase the solid 
surface density we performed computations of the Rosseland mean opacity at 1500 K 
for larger grain sizes. We increased the entire size distribution by a given factor y 
which means that we do not only grow large particles, but efficiently remove the small 
particles as well. Thus, where the MRN distribution we used so far runs from 0.005 fim 
to 0.25 fj.m, the new size distribution contains grains with sizes from (y x 0.005) /urn to 
(y X 0.25) yum. For these larger grain sizes the effects of anisotropic scattering become 
significant. The argumentation on the maximum total optical depth were made under 
the assumption of isotropic scattering. When anisotropic scattering is important, the 
same argumentation can be made when scaling the scattering cross section by a factor 
(I - g) (see Ishimaru 19781 where < ^ < 1 is the anisotropy parameter For values 
of g close to unity, most radiation is scattered in forward direction which for radiative 
transfer purposes is effectively a strong reduction of the effect of scattering. Values of 
g close to zero are more equivalent to isotropic scattering. 

The solid curve in Fig.|6]shows the increase in the maximum allowed solid surface 
density as a function of the scale factor y. For example, we see that by increasing the 
grain size distribution by a factor 1000 the maximum allowed surface density increases 
by a factor ~44. This is because at these grain sizes the Rosseland mean opacity (cor- 
rected for anisotropic scattering) is reduced by a factor of ~44. The dotted curve in 
Fig.|6]shows the same but now for very fluffy aggregates. The opacities here are com- 
puted using the Aggregate Polarizability Mixing Rule (APMR; see Min et al. 2008| l. 
It is clear that compared to the grain growth computed using compact grains (the solid 
curve) the fluffy aggregates have to grow to much larger sizes in order to reduce the 
opacity sufficiently to allow more solid mass to condense. If indeed sintering causes 
rigid bonds, prohibiting compaction of the formed aggregates, such fluffy aggregates 
are not an unlikely outcome of the aggregation process. 

An interesting feature for the compact grains (and to a lesser extend for the fluffy 
aggregates) in Fig.[6]is that the opacity of the grains first increases, which causes more 
mass to evaporate. Only when the particles are sufficiently large the opacity decreases 
and the solid surface density can increase. This is caused by the increasing scattering 
efficiency when the grains grow from interstellar grain sizes to several microns in size. 
It implies that in order to form large particles, a barrier has to be taken. The complex 
behavior of grain growth under these conditions asks for a more detailed study. For 
the fluffy aggregates this effect is much weaker. This is because for these aggregated 
structures the scattering is heavily forward peaked which makes the effect of scattering 
on the radiative transfer very small. 

Another possibility to decrease the opacity of the grains in the hot regions is to take 
away the species that contribute most to the opacity. We took a mixture of silicates, 
iron sulfide and carboneceous grains. It is well known that iron sulfide is only stable 
at temperatures below 680 K ( [Pollack et al.| [T994| . This means that the iron sulfide 
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Figure 6: The increase of tlie maximum allowed solid surface density as a function of the grain size distri- 
bution scaling factor y in the thermostat region of the disk. The solid curve is for the case of solid particle 
growth (computed using the Distribution of Hollow Spheres, DHS), the dotted curve is for the growth of 
fluffy aggregates (computed using the Aggregate Polarizability Mixing Rule, APMR). 

opacity source is most likely not there, though it can be argued that metallic iron will 
take over. Removing iron sulfide from the mixture we decrease the Rosseland mean 
opacity at 1500 K by a factor ~1.2, so this would allow only somewhat more mass to 
form. The dominant species of opacity in our computations is carbon. If we were to 
remove all the carboneceous species but leave the iron sulfide we reduce the Rosseland 
mean opacity at 1500 K by a factor ~3. Indeed, we find that the surface density in the 
inner region for our model without carboneceous solids is a factor of 3 higher that our 
standard model. Removing both the carboneceous species and the iron sulfide reduces 
the Rosseland mean opacity at 1500 K by a factor ~20. Though we did not perform the 
exact radiative transfer computations we can deduce that this would allow ~ 20 times 
more mass in silicate grains to condense in the hot inner regions of the disk. This is 
assuming no other species adding opacity, like for instance metallic iron, form in these 
parts of the disk. Since metallic iron is a very strong opacity source (even slightly 
stronger than iron sulfide) even a small amount would already increase the opacity 
again to approximately the same level as with iron sulfide and carbon present. 

In conclusion we find that the solid surface density in the region controlled by evap- 
oration due to viscous heating is completely determined by the opacity of the grains 
and thus heavily depends on the local grain size and chemical composition. 

43. Coagulation timescale in the thermostat region 

Another important effect of the opacity-temperature control system is that the re- 
duced number density of small dust grains will have a significant influence on the 
coagulation time scale in that region. The time scale at which coagulation can proceed 
depends strongly on the initial number density of grains at which the aggregation pro- 
cess has to start. Smaller particle densities increase the mean free path for particles and, 
consequently, the collision time for a dust particle to find another one to stick to and 



18 



grow. This eifect is extended through the further growth phases. Dominik and Dulle- 



mond ( 2008) have shown that this effect extends the time scales for removing small 



grains from the disk by coagulation approximately inversely proportional to the initial 
number density. This means that that coagulation will be strongly slowed down in this 
region, and that the time scales over which opacity dominates the number of grains that 
can be present in the thermostat region might be determined by the mass accretion time 
scale of the disk rather than the coagulation time scale. We believe that this effect can 
have significant influence on coagulation and disk evolution models. Further studies 
on this issue are necessary and will be the subject of a follow-up paper. 



5. Conclusions 



We successfully computed the solid density structure and snow line in the early 
Solar nebula for various parameters of the mass accretion rate and grain composition. 
This was done using full 3-D radiative transfer taking into account both the energy 
released by viscous heating as well as radiation from the Sun. The density structure of 
both the gas and the solids are computed self consistently using an a description for 
the viscosity and detailed evaporation physics for the dust and ice. We find that: 

• The location of the snow line is determined by the opacity of the grains and the 
mass accretion rate. We show that approximate radiative transfer methods to 
compute the location of the snow line are quite accurate. 

• At high mass accretion rates the midplane temperature in the inner few AU can 
exceed the silicate evaporation temperature. This causes a balance between con- 
densation and evaporation of matter acting like a thermostat keeping the mid- 
plane at the silicate evaporation temperature over a significant region. The solid 
mass in the inner regions is dramatically reduced by this effect and models com- 
puting the solid state surface density of the early Solar system should take this 
into account. 

• For our standard set of parameters the surface density in the solid state in the 
region of the terrestrial planets is below the MMSN for all values of the mass ac- 
cretion rate. This poses serious constraints on the scenarios of planet formation. 
It could point to a significantly lower opacity value at the time of planet forma- 
tion, which can be achieved by increasing the average grain size or significantly 
change the chemical composition. 

• At high accretion rates the increased midplane temperature of the disk due to 
viscous heating causes the midplane temperatures to be > lOOOK in a relatively 
large part of the disk. Under these conditions thermal processing of dust material 
as well as sintering of aggregates becomes important. 

• We find that for the high mass accretion rates the effects of convective cooling of 
the midplane are important to take into account. 
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